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Abstract 


The stability and Hopf bifurcation of a nonlinear mathematical model are 
described by the delay differential equation proposed by Wodarz for inter- 
action between uninfected tumor cells and infected tumor cells with the 
virus. By choosing 7 as a bifurcation parameter, we show that the Hopf 
bifurcation can occur for a critical value 7. Using the normal form theory 
and the center manifold theory, formulas are given to determine the sta- 
bility and the direction of bifurcation and other properties of bifurcating 
periodic solutions. Then, by changing the infection rate to two nonlinear 
infection rates, we investigate the stability and existence of a limit cycle for 
the appropriate value of 7, numerically. Lastly, we present some numerical 
simulations to justify our theoretical results. 
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1 Introduction 


Cancer is a significant cause of death in the world. Thus, it is essential to dis- 
cover some practical ways to prevail over it. Many studies have been made 
on cancer treatments, tumor cells behavior, clinical care, and so on. The 
primary purpose of cancer treatment is to reduce the destructive effects of 
cell behaviors [7]. The routine therapeutic substances for cancer are surgery, 
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radiation, and chemotherapies. The first two treatments are excellent choices 
for tumor cells when they have metastatic behavior [16]. The conventional 
therapies are not only efficient but also highly toxic, so most efforts focus 
on establishing tumor cells targeted treatments. For this reason, the biolog- 
ical encountering is the most up to date practice [11]. An Oncolytic virus 
is a type of virus that infects a cell and then explodes it. In this process, 
the cell dies, and new virus particles are produced. The main problem with 
Oncolytic virus therapy is that the virus stays in the blood for a short time 
due to the secretion of the immune system [6, 15]. Oncolytic viruses that 
specifically target tumor cells are promising anti-cancer therapeutic agents. 
Due to our lack of understanding of the Oncolytic virus’s dynamic spread 
in cancer cells, it is challenging to continuously control or eradicate cancer 
cells. The interaction between an Oncolytic virus and tumor cells, a type 
of virus-cell interaction, can be described by the mathematical models. The 
interplay between populations of uninfected tumor cells and infected tumor 
cells with the virus is complex and nonlinear. In this context, some mathe- 
matical models for the virus therapies of cancer cells can be seen as a tool for 
perceiving cancer-virus dynamics and finding better strategies for treatment; 
see [3, 4, 9]. 

In the proposed mathematical models, some parameters play a crucial role 
in the model’s qualitative analysis. For example, an average and optimal rate 
of virus-infected cell death may optimize the treatment success, or the lower 
the number of uninfected cancer cells in the stable state can better predict 
the treatment process; see [11, 12]. Here, we recall a set of mathematical 
models that describe the virus’s spread through the tumor cells in different 
ways [1]. First, we consider a general model as follows: 


d. 

< = 2F(2,y) — ByG(a,9), (1) 
d 

Th = ByG(a,y) — ay. 


This model consists of two populations: the population of uninfected tumor 
cells x and the infected tumor cells population by the virus y. The function 
F(a,y) denotes the growth rate of noninfected cells, and the function G(z, y) 
represents which the uninfected cells become infected with the virus. The 
coefficient @ represents the infectivity of the virus. Virus-infected cells die 
with a rate of ay. Assuming that the growth of tumor cells approaches the 
carrying capacity and after a while, slows down and that the cell growth rate 
reaches zero, F(x, y) = rx(1—(x+y)) can be expressed by the logistic func- 
tion. The G(z,y) function plays a crucial role in determining the system’s 
stability and helps us to make long-term predictions about treatment out- 
comes. Also, G(x, y) can be divided into two different classes. In the class [, 
if the number of noninfected cells is greater than the number of infected cells 
and the tumor cells are not solid, then the virus replicates and increases the 
number of infected cells. Biologically, the growth of the virus is exponential 
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and is called “fast virus spread.” In class H, when the number of infected 
cells increases, the virus’s growth rate decreases. The situation occurs in 
solid tumors because, as the number of infected cells increases, internal cells 
are surrounded by outer cells. Therefore, they cannot spread the virus. This 
type of infection is known as the “slow virus spread” [14]. 


We choose the following rates of infection given by [14]: 


Gi(2,y) =2, 

.. (bee 
Go(z,y) = (a@tyte)’ 
G3(z,y) = Gike 


Here G and G2 belong to class of “fast virus spread” and G3 belongs to class 
of “slow virus spread”. By replacing F(z,y) = rz(1 — =") and G = G, in 
(1), we get 


dx x+y 


= ra(1— )— Bay, (2) 
d 
7 = Pry — ay. 


In 2016, the following model was developed by Wodarz, similar to (2) in 
which tumor cells death rate was considered (see [13]) as 


a ra(1— ) — da — Bay, (3) 
dy 
Ge = Pty — oY 


so that 6 is defined the death rate for the uninfected tumor cells. Since (3) 
is an extension of (2), we introduce a delay time 7 for entering virus in (3), 
such that it changes into 


) — da — Bay, (4) 


Because of the simplicity of the delay calculations in system (4), first, we 
perform the calculations the delay equation of (4) analytically, and we find 
the appropriate parameter of 7 to produce the limit cycle in (4). Now, by 
changing the linear infection rate of (4) to Gz and G3 and entering 7, we 
obtain, respectively, 
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a all +e)ey 
a =ra(l1 ; )— da ere? (5) 
dy  B(1+e)a(t—r)y(t—T) ‘ 
dt (a(t-—7r) + y(t —T) +6) mt 

and 
dix tty, s. Bay 
dt 7 oe W ° (xy3 oe)! (6) 
dy _ __ Ba(t—T)y(t—7) a 


dt (x(t —r)y(t—7)3 +e) 


By placing the parameter 7 obtained from (4) and changing €, we simulate 
the existence of a limit cycle in (5) and (6), numerically. The parameters 
r, w, €, 8, and a are all positive and independent of the time. The state 
variables x(t) and y(t) are nonnegative. The parameter r is the reproduction 
rate of tumor cells, w is the maximum carrying capacity of tumor cells, and 
the coefficient 3 denotes the virus’s replication rate. Moreover, a is a death 
rate for the population y that are killed by the virus and €¢ is a positive 
parameter that is expressed only to fit the data to the model better. 
In [10], the mathematical model 


dx ut+y 
x 


a — —b 
cree x) ~ Oty, 

d 

7 = ba(t— r)y(t—T)e-”” — ay, 


for tumor virotherapy with the viral life-cycle is considered. In this model, 
the delay parameter displays the period of the viral life-cycle. The two pa- 
rameters, b and a, are dominant factors in virotherapy. When b < a, the 
tumor cells reach their maximum size and the equilibrium solution (1,0), for 
T > 0, is stable. So, the therapy fails, but for b > a, the model’s dynamic is 
much complicated, where the viral life-cycle comes to play an important role. 
As b> a, the viruses cannot exterminate the tumor cells without the viral 
life-cycle period. When the delay parameter value is small, the tumor cells 
and the viruses coexist, and their equilibrium solution is locally asymptotical 
stable. If the value of the delay parameter is longer than the period of viral 
life-cycle, the coexisting solution will be unstable, then the population of the 
tumor cells and the viruses will not rest on a fixed level. The model has 
a stable period solution for the value of the delay parameter in the middle 
range. 

In appearance, by supposing 6 = 0, the model of (4) is similar to [10] 
for n = 0, but they are different in some ways. Firstly, in (4), there is the 
mortality rate of uninfected cells, which results in different reproduction rates 
with [10]. In this article, we investigate the stability of equilibrium points 
of (3) for different values of the reproduction rate by using the Lyapunov 


Hopf bifurcation analysis in a delayed model of tumor therapy ... 163 


function. Then, with the Dulak theorem, we prove the absence of a limit 
cycle in (4) for 7 = 0. Another difference between this work and [10] is 
the use of the center manifold theory in the study of the stability, direction, 
and period of periodic solutions of bifurcation at (4) for critical values of 7, 
analytically. 

In [2], another mathematical model 


dx x+y 

ge el )- dx — Bry, 

dy x+y 

ae 7 Pty + sult 7 )—ay(t—7), 


for tumor virotherapy is proposed. Wodarz in his article has stated that the 
infection rate and death rate of virus-infected cells play a key role in model 
dynamics. In that article, there is the time delay in death of infected cells 
and it entered in the linear sentence. Also, the Hopf bifurcation for linear 
infection rate is investigated. Our article’s innovation is that the time delay in 
the arrival of the virus to tumor cells, and it entered in the nonlinear sentence. 
Also we have used three models with different Wodarz infection rates, two 
of them belong to the class of nonsolid tumors, and the other belongs to 
solid tumors. Using the obtained value 7 suitable for the existence of a limit 
cycle in the initial model and its placement in the next two models, we have 
numerically shown that this dynamic depicts the general biological conditions 
of solid and nonsolid tumor cells and the effect of viral therapy on them. 

In this article, we study the stability and Hopf bifurcation of system (4); 
then we compare it with (5) and (6). At first, in Theorem 1, we investigate 
the stability of the equilibrium points for system (4) for 7 = 0. Then in 
Lemma 2 and Theorem 3, we consider the existence of periodic solutions and 
the Hopf bifurcation for system (4) by inserting 7 into the second equation (3) 
as a delayed time for importing Oncolytic viruses to the body and countering 
it with cancer cells. In Section 4, we describe the stability, direction, and 
period of the bifurcating periodic solutions at critical values of 7, by using 
the center manifold theory introduced in [5]. In Section 5, we substitute the 
appropriate value of 7 obtained from (4) into (5) and (6). Then we simulate 
(5) and (6) for different values of ¢. Finally, we discuss about them. 


2 Stability of equilibrium points 


In this section, we obtain equilibrium points of system (4) for 7 = 0 and study 
the conditions for the existence of positive equilibrium points. It is obvious 


that system (4) has equilibrium points Eo = (0,0), Ey = (29,0), and 


—_ fa Bw(r—d)-ra 
By = (4, B(r+Bw) ). 


The Jacobian matrix of system (4) is denoted by 
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low) =(" 5) s (4+ B)y ae, 


Theorem 1. Suppose Ro = 5; then for 7 = 0, we have the following prop- 
erties: 


(a) If Ro < 1, then Eo is a unique equilibrium point that is a stable node. 


(b) If 1 — = < 3, and Ro > 1, then there exist two equilibrium points Eo 
and FE, such that Eo is a saddle point and FE is a stable focus. 


(c) If 1— Ro > gy and Ro > 1, then there exist three equilibrium points 
Eo, Ey, and E2 such that Ep and FE, are saddle points and EF is a stable 
focus. 

(d) If 1— a = a then E£, = E> and we have two equilibrium points such 
that Eo is a saddle point and £| is locally asymptotically stable. 


Proof. 


(a) The Jacobian matrix of system (4) at the origin is 


J(0,0) = Po (7) 


which has the eigenvalues \; = r—d < Oand A; = —a < 0, when Ro = . <i. 
This shows that Ep is a stable node. 


(b) When Ro > 1, for 7 = 0, there exists the equilibrium point £, and the 
Jacobian matrix of system (4) at FE is 


r4 r p ee 
J(E1) = ( 0 ? ( pean m ) : (8) 


Tr 


It has the characteristic polynomial 
M+ mA+ Me = 0, 


where 


For a > Bw Bus or 1 = < Ba we have that my, > 0 (k = 1,2) and this 
implies that FE, is a stable focus. 

(c) For Ro > l anda < Bw— Bu, there exist two positive equilibrium points 
FE, and E»2 such that the Jacobian matrix of system (4) at E> is 
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are) 
J (Es) = = . : 


r+Bw 


Since tr(J(E2)) = Zot < 0 and det(J(E2)) = oot —ra) > 0, then Ep is 
a stable focus. 


(d) By substituting 1 — Re = a in (8), we obtain 4; = —(r — 6) < 0 and 
A2 = 0. Hence, FE; is a nonhyperbolic fixed point with a zero eigenvalue. 
To show its stability, we use a suitable Lyapunove function. In the set A = 
{(x(t), y(t))|x(t) > 0, y(t) > 0}, we define a Lyapunov function 


w(r—d) w(r—6d) m2 _ Bw tr 


r r wir—d)' Bu a 


V(a,y)=2 


such that V(£,) = 0 and V(z,y) > 0 for all (x,y) € A\ {Ei}. We now 
show 4 < 0 for all (x,y) € A. By differentiating V(x, y) along the solutions 
(x(t), y(t)) of (4), we have 


wa ES+ Fe -( wen) 4 Ey 


dt Ox dt dy dt Bu 


(: ura) G ae — Ta? yeu) ! (EE) ay ay) 


Ww 


= 2(r — 6) x + (r+ Bu) (=) i) 


w [r? r r+ Bw Bu(r—6)—ra 
a 27(r d)at+(r ay?) +( Bus )( ( 2) )y 


_ =a (me (r 9) : (5) (eee), 


Asl—g = Zt, then #9" _ 0 and = # (fe (r 5)? <0. 


r Ww 


If we set on =0, then x = caleal and by substituting in system (4), we get 
y =0. Therefore, EF; is asymptotically stable. 


Theorem 2. System (3) in the set A = {(a(t), y(t))| x(t) > 0, y(t) > 0} 
does not have any periodic solutions. 


Proof. We define F, = (r — 0)x 


x (THB#) oy and Fh = Bry — ay. By 


Tr 
Ww 


considering the Dulac function L(a,y) = re we have 
O(LF;) O(LF 2) = r <0. 
Ox Oy wy 


Thus, by the Bendixson’s criterion, system (3) has no closed orbits lying 
entirely in A. 


We have now shown that system (3) has no periodic solutions. We con- 
sider system (4) and show that the Hopf bifurcation exists under special 
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conditions for tT > 0. First, we suppose that E* = (x*,y*) is a positive 
equilibrium point for the system (4). The characteristic equation of any 
equilibrium point E* = (x*,y*) of system (4) for 7 > 0, is given by 


D(A,T) = det(AI — Qi — Qe) = 0, 


where 


_ AN-—r +6474 (£4 B)y* (24+ 8)z* 
D(A,r) = aet ( =pyre-** A+ pre" La 


=)? } A,X } Ao (A { B,)Boe~", 


where 
7) i ae x 
Ay =a (r 6) By ’ 
Ww Ww 
2ra* ry” 
Ag = 6) 4 + By* 
p=a((r—6)+—— + 4 By"), 
2 * 
Hae 
Ww 
Bo = Bxu* +a. 


Lemma 1. Let x«(0) = ¢1(0) > 0 and let y(@) = ¢2(@) > 0, for @ € [—7,0], 
where ¢@ = (¢1, ¢2) € C([—7, 0]; R?). Then the solution (x(t), y(t)) of system 
(4), defined on the interval A = [0,7] for some 0 < T < ~, is positive and 
uniformly bounded on A. 


Proof. The first equation of system (4) is a Bernouli equation, which can be 
written as 


a(t) — a(t) {(r — 6) - (+ B)y(t)} = -22°(0. 


Ww 
Its solution with the initial data x(0) = Zo is 


-1 


t 
a(t) = pedo Lr-9) (E+ By ay {0 f F elo {(r-I) (Et AVM far gg 4 i} 
0 W 


Thus, x(t) > 0 when z(0) = 2p > 0. We know y(6) > 0 for 6 € [—7, 0]; then 
y(t — r) > 0 for t € [0,7]. Thus, by the second equation of system (4), we 
have y(t) > y(0)e~™ > 0, for t € [0,7]. In a similar way, we obtain y(t) > 0 
for t € [7,27]. By repeating this process, we conclude that y(t) is nonnegative 
for all t > 0. To prove the boundedness property of the solutions, we will 
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first show that x(t) is bounded. By the first equation of system (4), we have 


#(t) < (r — d)a(t) — “a(0), 


and hence P 
a(t) < a =o 
pa: (es _ r) en (rat 
Putting 8) = a —r, we have 
a(t) < w(r — 0) 


r+ BoeW Ht , 
and therefore lim sup,_,,, x(t) < wr 9) Now by supposing Q(t) = x(t) + 
y(t +7) and y =r—06 > 0, we have limsup,_,,, x(t) < “* =: 6; and 


= yx(t) — —2?(t) — —a(¢)y(t) — Bx(t)y(t) + Ba(e)y(t) — ay(t + 7) 
= y(t) — <27(t) — <a(t)y(t) — ay(t +7) — a(t) + ax(t) 

< f(y +a) ~ © 0?(t) — “a(t)y(t) — ay(t + 7) — ax(t) 

< Bly +a) — aly(¢+r) + a(t) 

< Bly +a) — aQ(e). 


Hence, Q(t) + aQ(t) < B1(y + a), which gives Q(t) < enka) (Ante) 


Q(0))e~*. This implies that limsup,_,., 2(t) < eu) forr > 6. Asa 
consequence, the functions x(t) and y(t) are uniformly bounded. 


3 Hopf bifurcation 


We claim that for some 7 > 0, a Hopf bifurcation occurs. The interior 
equilibrium point E> exists when Bw(r —6)—ra > 0. We move this point to 
the point (1,1) by setting 


L= ma 
Bw(r — 6) — ra 
2 B(r + Ba) vu 
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Applying the transformation (9) to (4), we get 
dX _ 


aT = X(y —mX) — nXY, (10) 
dy 
aT = aY +aX(T—7)Y(T—7), 
where Bur — 5) 
ra w(r—d)— ra 
yer-6 maT ne Ea (11) 


Therefore, (X,Y) = (1,1) is an equilibrium point for system (10). Writing 
t instead of T and considering ui(t) = X(T) — 1 and we(t) = Y(T) —1 in 
system (10), we find that 


ty (t) = (ui (t) + 1)(y — miu (dt) +: 1) — n(ua(t) + 1)), (12) 
u(t) = —a(ue(t) +1) + a(ui(t — rT) + 1)(we(t — 7) +1). 


From (11), we know y—m-—n = 0, and by simplifying the right side of (12), 
the linear part of system (12) around u = 0 is given by 
u(t) = —mur(t) — nua(t), (13) 
ua(t) = —aue(t) + aui(t — rT) + aua(t — 7). 


We know that the characteristic equation for system (13) is given by 


D(A, T) = det(AI — Qi — Qoe~) = 0, 


—m —n 00 
a=(4 mF a= (08): 


By the above assumptions, we can write 


where 


D(A,T) =X +pirA + p2- (A+ a)ne’, (14) 


in which 
Pi=M+a, pz =ma, ga=m—n, qg =a. 


For tT > 0, we show that under some conditions on parameters, a Hopf 
bifurcation happens. We suppose that 7€ for some € > 0 is a root of (14) so 
that 
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D(iE, 7) = (i€)? + (t€)p1 + po — (1€ + 1) ne" 


—€? + ip, + po — i€qe cos(Er) 
— qig2 cos(ET) — Eq2 sin(€T) + igige sin(ET) 
= 0. 


Separating the imaginary and real parts of (15) gives that 


pz — €? = qiqe cos(Er) + Eqz sin(Er), 


Epi = €q2 cos(€r) — qigz sin(ér). 
Eliminating sin(é7) and cos(€r) from (16) and (17) implies that 
C+AC+B=0, 


where 
A=pi-2p-@, B=p3- aid. 
By assuming Z = €7, equation (18) changes into 
ZsrAg ea =|, 


which has the solutions 


£7 + i€py + po — (i€ + G1) q2(cos(E7) — isin(E7)) 


—(pt — 2p2 — g3) + V/(p? — 2p2 — G3)? — 4(p3 — G3) 
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(15) 


(16) 


(17) 


(18) 


(19) 


442 = 5 
We suppose that h(Z) = Z? + AZ + B has two positive roots Z,, k = 0,1 
such that 
0< 4 < %, 
and 


fo =V2Zo, 


Then we have the following result. 


&=J/Z. 


Lemma 2. Suppose that h’(Z,) #0 and Z, = €. Then dr (0) # 0 and 


its sign is given by the sign of h’(Z;,). In particular, we have 


h' (Zo) <0, h'(Z1) > 0, 


such that 


dr(7”) £0 dr(r\??) 
dtr , dt 


>0, f=0,1,2,.... 


Proof. Multiplying (16) by q; and (17) by €, we have 


(20) 
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pon — €q1 = qq cos(Er) + €qide sin(Er), 
pie? = qn cos(Er) — Eqiqz sin(Er), 
and hence 


peg — €2q1 + pié? _ peg — 67qi + pik? 
Gg + 6G qo(qi + €?) 


cos(€T) = 
Consequently, 


2 2 

(i) _ 1 | (2 =a +e )| , 

T,’ = — |arccos +2jn, for jeEZ. 
B  & qo(q7 + €?) 


Therefore, tig, for k = 0,1 are two pairs of imaginary eigenvalues for T = 
70). 
Let A(r) = r(r) + i€(7) be a root of (14). Then for rT = 7, we have 
Mr) = r(7) + ig(7) such that r(7) =0 (7) =&>0 
kk) = TT, k kt) = 9%; ko) = Sk > U. 

We substitute A(T) into (14), and we get 
d? (7) + pi A(T) + po — (A(T) + 1) qe” = 0. 


Differentiating this equality with respect to 7 leads to the identity 


(2X(7) 4 po (1 r(Mr) + 1) pet BOD 


T dt 


+ Ar)(A(T) + qr)qze*”” = 0. 


This gives 
(2a) " _ L 2A(r) + D1 7 
dr MAT +H) ATA +q)QerO™ — -A(7) 
X(T) 2X? (7) + p,A(T) T 


M(TV(A(T) Fa) M(T)(A(T) + a) qae AO) ACT)" 


Now, by using (14), we obtain that 


eas _ —U (7) — pa r 
dr M(T)\(A(T) Fan) AP(T)(A2(7) + pi A(T) + p2) A(T)’ 


and therefore, 
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(etn 


rar) 
- (Re sasat + a| — Percaesene val 
Rel) Ley 


Since, Nr) = r(7) + i€(7) = 1&;, hence 


—H _ q(—ikk +) 
“ (eat ss a) a lace Hale + a| 
2 
EER +i) 
p2 — A?(r) ) | p+ | 
a (emorn + piXA(T) + pa) J Irae me —€2(—@ + i€kpi + p2) 
(p35 — ER) 


En((—Gj + pa)? tise) 


—T 
rar) Fa 


By the above relations, we get 


een 


Let & = Z,. Then 


ay 


Since Z;, is a root of (19), we have 


_ 4 Pa — &% 
rar GREP +g?) E2((—€2 + po)? + E2p?) 


a ps — Zi, 


rat? Ze(Ze+ ai) Ze((—Ze + pr)? + Zevi) 


(—Z;, + pz)? — Zep, = Ze + (pi — 2p2)Ze +75 = BZe + GG = G(Ze+G), 
Ze + (pt — 2po — 43) Ze = —p3 + G7 03- 


Consequently, 


nen 


Ze-Pa+ GG — 2Zp + (pi — 2pe — 43) Ze 


r=1 — Zq3 (Ze +7) Zn93 (Ze + G) 
_ 22, + Pi — 2p2-@B _ h'(Ze) 
g3 (Ze + G2) Ty 


where 
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T; = g3 (Ze ale qi) = 0, 
hi! (Zp) = 2Z,, + p} — 2p2 — 5. 


Therefore, 
dr(r)]~" dReNr)] > | 3) ee 
sign | = I. 7 sign | an vent? = sign T, | > sign [h' (Zz) 


Theorem 3. Let (19) have at least a positive root. Then system (4) has a 
Hopf bifurcation at 7 = rf ) and has one periodic solution surrounding E*, 


where : 7 
7) = — fos (Be 5 a = )| + 29m. 
V Z* qo (qq +Z ) 


Proof. This result easily follows from Lemma 2. 


We recall the results stated in [8] on the roots of (14) for 7 > 0. 

Lemma 3. ((8]) Consider the following assumptions: 

(Al) pi — @ > 0; 

(A2) p2 — G92 > 0; 

(A3) pi-2p2—-q3>0 and py—qiqy>O0 or (pj — 2p, — 3)” < 
A(p3 — 97.43); 

(A4) pi > 2p2— 92 <0 or pa—agigg <0 and (pj — 2p2 — 3)” = 
4(p2 — 9199); 

(A5) pi — 2p 92 <0, pe — aig > 0 or (pi — 2pa — a3)” > 
4(p3 — 4192). 

Then the following properties hold: 

(i) If (A1) — (A3) are satisfied, then the real parts of all roots of (14) are 


negative for T > 0. 


(ii) If (A1), (A2), and (A4) are satisfied, then (14) at 7 = 7) has a pair of 
imaginary roots +7, such that the real parts of all roots except +i&; 
are negative. 


(iii) If (A1), (A2), and (A5) are satisfied, then (14) at 7 = 7), has a pair of 
imaginary roots +2&;, such that the real parts of all roots except +7£; 
are negative. 


Theorem 4. If there are no positive roots for h(Z) = Z*? + AZ + B and 
Ro > 1, then Ey is locally asymptotically stable for any 7 > 0. 


Proof. Lemma 3 states that the real parts all of roots of (14) are negative, 
and this shows that FE» is locally asymptotically stable for any 7 > 0. 
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4 Stability of the Hopf bifurcation 


In the previous section, we derived some conditions on Hopf bifurcation at the 
equilibrium point E* = (1,1), when 7 = 7, (j = 0,1,2...). In this section, 
we obtain stability, direction, and period of periodic solution bifurcating from 
certain values of tT. We use the center manifold and normal form theory 
recognized in [5]. Let system (4) have a pair of imaginary roots when tT = T* 
and let the system have a Hopf bifurcation from E*. We set 


ui(t) = X(rt)—-1, ug(t) = Y(rt)—-1, THT +y, VER. 
Then, from (10), we get 


ta(t) =(7* +) (uf) + 1)(—mui (4) — nua(t)), (21) 
ug(t) = a(r* + v)(—ue(t) + u(t — 1ue(t — 1) + u(t — 1) + ue(t — 1)). 


Equation (21) can be written in the form 
u(t) = Ly (ue) + f(v, ut), u(t) = (ui(t), u2(t)) € R’. 


Let C = C({-1, 0]; R?) be the Banach space of continuous functions on [—1, 0] 
with values in R?. Then for ¢ = (¢1, 62) € C, we define 


Lyx PR, f:RxC—R?’, 


through 


ronieon( 2) 
flv, 6) = (7 +) [=m Eo -_ 


_ (n* —m¢7(0) — n@1(0)¢2(0) 
= C +) ( i, (22) 


By the Rise representation theorem, there exists a bounded variation func- 
tion 7(0,v) for 6 € [—1, 0], such that 


0 
n= i dn(6,v) (6), 


= 
with 


m —n 


n(0,v) = (r* +0) & ) 5(0) — (7* +) (° 0) 5(0 +1), 


where 6 is Dirac delta function. For ¢ € C'([-1, 0]; IR), we define 
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dp(0) _ 
ene {0 dn( (6, s)(s), : - : (23) 
and ; 7 
ee ty f(v.9),9 ay : (24) 
By the above notations, system (21) can be written as 
ue = A(v)ur + R(v)ut, (25) 


in which u;(0) = u(t + @) for @ € [-1,0]. For ~ € C1((0, 1]; (R?)*), define 


dep(s) 


arone)={ haa se (26) 


(t,0)¢(—t), 8 


and denote the bilinear inner product by 


0 0 
< (8), 0(8) >= 0(0)9(0) — I. - BE — A)dn(A)o(E)dg. (27) 


We suppose that 7(@) = 7(0,0) and that A* is the adjoint operator of the 
linear map A(0). From Section 3, we have that +i€*7* are the eigenvalues of 
A(0) and that also the eigenvalues of A*. Assume that q(@) is an eigenvector 
of A(0) for i€*7* and that q*(s) is an eigenvector of A* for —ig*r*. Then by 
the definition of eigenvalue and using (23), we have 


A(0)q(@) = i€*7*q(0) = d (8) =i€*r*q(0) for  6@€[-1,0), 


and hence 
q(8) = q(0)e 7", = -1 5 0<0, 
q(-1) = 40 ser 
Aoq(0) + Bog(—1) = 1€*1*q(0) 
Thus, 


(Gree —a— ito tines) (at) 7 & , 2) 


It is easily seen that a nontrivial solution of (28) is as follows: 


350), (0) =140, 


so that 


g(8) = (41(0), 42(0)) "7° = (1. ae) “ied, 


n 
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Similarly, we obtain q*(s) = E (q%(0),q%(0))* e®°7"S as an eigenvector of A* 
corresponding to the eigenvalue —ié*7*. Indeed, by using (26), we have 


A*q*(s) = —i€*7*q"(s) = —q'i(s) =—i€*r*q(s) for 
and therefore 
q*(s)=q*(O)e*7™*, O<s<1, 


q*(1) =q*(O)e*™, 
Agog" (0) + Bog" (1) = —i€*6"q" (0). 


Hence 


Equation (29) shows that 


. sci 
if* -—at+ae’™ , 


qi (0) = q3(0), qi (0) =140. 


n 


Thus, we find that 


€ (0, 1], 


sect ig* — gtr” Sank 
a"(8) = B (af(0), a5(0)) "7 ‘=2(* << 1) 9 


n 


For normalizing g and q*, we assume that < q*(s),q(0) >= 1 and by using 


(27), we get 


< q°(s),q(9) >=E(F*1(0) he q2(0))" 


-[, f= El 1(0), Le Fd) (1, 2(0) dé 


ee ee —m — i€* 
a ee eee Gee 


0 0 ek —it* 7* ae 
-| | ( ig* — a + ae 2) eT (E-9) dn (Q) 
—1/€=0 th 


¥ (1 ss ) eft" 8aet 
nm 


=H{ —2i€* —a+ae—*"™ —m 
n 


0 _ gee —ig*r* ‘se. 
/ ( Sa a) des 7° dn(8) ( is 
Ls n os 
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=E 


ro eae Tm re Arn) Ha" aero} 
n n 
=1 


d 


which gives 


E= 


= {—26" —a+ae—®'T —m—r*(a(m—n) + taé*)e y 
i 


We now describe the solutions of (21) on the center manifold Cp at v = 0. 
Let u(t) be the solution of (21) at v = 0. Then 


u(t) = uo(t) + U(uo(t)), 
for uo(t) € Eo and Eo = span{q(6), q*(6)}. Thus 
uo(t) = 2(t)q(9) + 2(t)q" (8) = 2(t)q() + Z(t) a), 
and hence u;(#) on the center manifold can be written as 


ux(9) = z(t)q(@) + 2(t)q(8) + U(z(E), 2(¢)) 
= 2Re{z(t)q(@)} + W(t, 8), (30) 


where 

W(t, 0) = ur(@) — 2Re{z(t)q(A)}. 
It is clear that if u:(0) is real, then W is also real. Now, we set 
2 3 


Bs | (areal ae (31) 


5 a 
W(t, 0) = W(z, z,9) = Wao + W122 + Woo 6 


in which z, Z are the coordinates of uz with respect to q,q*. The reduced 
equations up to order three for (z,Z) will be computed in the following. 
Equation (25) used with v = 6 = 0 yields 


u(t) = A(0)u(t) + R(O)u(t), 
which is, by (30), equivalent to 
3(t)q(0) + Z(t)q(0) + W(t, 0) = A(0)u(t) + R(0)u(t). (32) 


Since the inner product of the left-hand side of (32) with q*(0) is equal to 
2(t), then 
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2(t) =< q*(0), A(O)u(t) + R(O)u(t) > 
=< q*(0), A(O)u(t) > + < g*(0), R(O)u(t) > 
=< A*(0)q* (0), u(t) > + < q*(0), R(O)u(t) > 


where 


2 2 222 
g(2(t), 2(t)) = 920-5 + 9112Z + go2 5 + gai fore, (33) 


To simplify the computations, we let 


—m — i&* 5 i€* —a + aes’ 
n0)=—"—* =u, gy =" =N. 
n n 
Then (30) reads as 
ue(8) = 2(t)q(9) + 2(t)q(8) + WE, @) 
= (1,M)Te®"9 a(t) + (1, M) Te "7 z(t) 
2 2 
+ Wao 5 +t Wi12z + Wo2 5 ae (34) 


It follows from (34) that 


@)=24+274+W™ (0), 

(t) = Mz+Mz+Ww”(0), 

u(t —1) = ze" 44 ze*®T™ 4 WY (-1), 
ug(t — 1) = Mze*"" + Mze*™" 4+ W(-1), 


U1 
U2 
t 22 


so that W(t,0) = (W (6), W@)(6)) € R?. By replacing ¢(0) = u,(@) into 
(22) for v = 0 and the above relations, we have 
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9(%,2) 
=" (0) f(0, ut) 
aa ( —m3 (t) — ns (Qu(0) 
=q"(0)r ( ae aa nade, ) 


—m(z + 2+W)(0))? — n(z +2 + WO (0))(Mz + Mz + W®)(0)) 


_F(N,1)r* 
a(ze— i867" | zeié*7™ + W()( 1))(Mzei6"7" + Mzets*t* + W)(-1)) 


=— ENr*{m(z+2+W) (0))? + n(z+ 2+ W (0))(Mz + Mz + W)(0))} 

+ Er* {a(ze78 7" + ze8" 7" 4 w( 1))(Mze7#8"™" + Mze's**" + W)(—-1))} 
52 

z 


2 
=~ ENr*{2(m+nM) > + (2m + nRe{M}) 22 +2(m+nM) 5 


42 (2mw{? (0) + 2mw) (0) + nMW? (0) + nw) (0) + nw? (0) 
=2 
z 


2z oe 2 eee 
} nW39) (0) =" } + Br*{2 (ame ?"7") = + (aRe{M}) 22 +2 (aMe™"7") 


4. 2(ae#€°"* WAP (-1) + ae®"7" WP) (-1) + aMe—**"7" WA) (-1) 


4: aMfe*"*" w39’(-1)) | fi (35) 
We now replace (33) in the left-hand side of (35), and we get 

g20 = —-2ENr*(m+nM) + 2Er*(aMe—2*"7"), 

gu = EN7r*(2m+nM +nM) + Er*(aM +aM), 

goz2 = —-2ENr*(m+nM) + 2Er*(aMe?*' 7"), 

go = -2ENr*{(2m + nM)W4? (0) + (2m + nM)W4q) (0) + nW4y? (0) 
+ nWSD (0)} + 2Er*{ae~®°" WY (-1) + ae®™ WY (-1) 
+aMe—* wY (-1)+ aMe*’*” ws) (-1)}. (36) 


As Wi1(0) and W20(@) are unknown, we should calculate them. To this end, 


we consider 
W = t%& — qz — G2, (37) 


with 
tie = A(O)u, + R(O)u; = (A(O) + R(O))(qz + G2+W), 
2=1€*7*24+ GF (0)fo, 
B= iz + 4° (O)fo. 


Then (37) becomes 
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W =A(0)(W + gz + G2) — (i€* "2 + (0) fo)a — (—1€*7*Z + 4° (0) fo) a 
+ R(O)ut 
=A(0)W + A(0)qz + A(0)Gz — (i€*7* z)q + (1E* 77 ZG 
— 7 (0) fog — 4° (0) fod + R(O)ue 
=A(O)W + (16" 7" q)z — 1G" 7" Gz — (18* 7" 2) + (1b *7*2)G 
— 7 (0) fog — 4° (0) fod + R(O)ue 
=A(0)W — ¢°(0) fog + G* (0) fog + R(0)ue 
=A(0)W — 2Re{q (0) fog} + R(O)ue. 


(38) 


Next, from (24) and (38), we find that 
AW — 2Re{q (0) fog} + F(0, ut) 0=0, 
= AW + H(z, Z,9), (39) 
where 
2 = 
H(z, 2,6) = Hoo(8) > + Hi (0)2z2 + Ho2(8)> + °° : (40) 


By replacing (31) and (40) into (39), we obtain that 


Wao (0) 22 + Wi1(0)2Z + Wi (0) zz + Wo2(0)zZz sass 


bo 


2 z2 
= A(Wa0(8) = + Wi1(0)22 + Wo2(4) >) + Hoo (0) > 
z2 
+ Ay (0)22 + Hoo (0) peer, 
Also, 


Wa0(0)z(t€*7* z + g(z, Z)) + Wir (0) Z(i€*7* z + G(z, Z)) 


+ Wi1(0)2(—t€" 7" 2 + g(z, Z)) + Wo0(0)2Z(—1E* 7" Z + g(z, Z)) 
= (AW 9 (0) + Hy9(0)) = + (AWi1(6) + Ay,(0))2z 
+ (AWo2(60) + Hild) Se 


iw) 


Thus, 


2 ge 
(21 7* Wao (8) > + (—216*7* Woo (8) > +> 
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Comparing the similar terms in both sides of equality (41) leads to 


AW 20 (9) + H20(8) = 21" 7" Wa9(8) => (A — 216" 7")W20(8) = —H20(9), 
(42) 
AW1(0) + Ay, (0) =0 = AWi1(6) = —Hy,(6). (43) 
From (39) for 0 € [—1,0), we have that 
H(z, 2,0) = —2Re{q"(0) foa(9)} 
= —F' (0) foa(@) — (0) foal) 
= —9(2, 2)q(9) — g(2, 2)q(8). (44) 
Then by comparing the coefficients of (33) and (40), we get 
H9(9) = —g209() — Go24(9), (45) 
Ay, (9) = —9114() — G11 (9). (46) 
By substituting (45) into (42), we obtain 
AW 0 (8) = 21€* T* Wo0(0) = H49(0) 
= 21€"7"Wa0(8) + g2oq(4) + go29(4). (47) 


By the definition of A(v) for @ € [—1,0), we have AW29(0) = Wao(0), and by 
replacing it into (47), we get 


Woo (0) = 2%€*7* W9(0) + g209q(0) + Go2G(9) 
= 2é*7*Wo9(0) + goog(O)e*'*? + Go2q(Oe*™ 8, 


where q(@) = q(0)e**7 °. Hence, 
ok Ok a ok ok ook oe ok 
Wa9(8) = e°7"? 1 (s200(0)e'* 78 4 Goog(O)e~** 7 eee re ig + a 
) 
_ 1920 igtr*9 , 2902-1 —ig*r*0 2ie*r* 0 4 
&*r* q(O)e + 3E*7r* q(O)e ar Cie : ( 8) 


Similarly, it follows from (43), (46) and the definition of A(v) that 


AW11 (9) = —Hi1 (9), 
Wir (8) = grrq(0) + G119(9) = griq(O)e 7? + Gug(O)e" 7”. 


Then 
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6 
Wiu(8) = tf (guia(O)e*" 9 4 G11 G(O)e* “ dé + ct 
0 


tigi ie*r*9 , WIL 7 —ietr*O 
=— 0 —— q(0 Co. 49 
1881 ge"? + HO Goer"? 4 05 (49) 


We now need to find appropriate C, and C2 for (48) and (49). From the 
definition of A(v) for @ = 0 and (42) we have that 


0 
/ dn(@)W29 (8) _ 2iE*T* Wao (0) _ H9(0), (50) 
-1 

where 7(0,0) = 7(6). From the definition of H(z,Z,0) at (39) and (44) for 


0? = 0, we have 


29 (0) = —G* (0) foa(0) — 4° (0) fog(0) + F(0, ue) 


2(m+nM) ) 


= —g204q(0) _ 9029(0) te (ape (51) 


Now, we substitute (48) and (51) into (50), and we get 


r re (2 g(0) + 3 a0) + C1) 


—a Cor 3E*T* 
ef 2.0 7.920 —ietr* , 2902 _/\ sete" —2ig*7* 
T (", _) (= q(0)e _ 3e87* q(O)e +Cye 


= 2ig*r* (EP 9(0) + s-a(0) + C1) + 920410) + G00410) 


+rt( ae ). 


—2aMe-216°" 


It reduces to 
« ¢ 4920 —m —n 7902 —m —n _ 
EP (ee a pacers + Oe eee gy eee ) 10 


—m™m —n 
+ ( gear a+ aen 2 ) Cy} 
2(m +nM) ) 


= —g20q(0) — Go2G(0) + 2%€*7"C, + 7° (Soe 


Then we use the characteristic (AJ — Q1 — Qz2e~*")q(0) = 0 and obtain 
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* 1920 + ok 7Go2 ek = —m™m —n 
T {= (i€"q(0)) + 3e+ 7 i€*q(0)) + ( geeer a+ ae" ) a} 


_ oe ee J 2(m+nM) 
= —9209(0) — Go2q(0) + 21G*7"*C1 + 7 fo) 


Therefore, 


e - 2 M 
( a Des ) Ch = pu 7C, + ( agen ) 


ae—**"T" _a94 ae —2aMe~2is"7" 


which is equivalent to 


ede =n oY \ _ ( 2m+nM) 
aew26 7" q+ qe72tFT* — 2ig* of © \=2aMe*S" ]° 
It follows that 
CM = —2(ame7?"8"" — (a + 2i€)(m + nM)) 
1 ae~?*§*7* (ME +m —n) — WE(a + 2E+m) —am’ 


2a(m+nM — 2Migé—mM) 
a(2ié + m —n) — (2iE(a + 2E + m) + am)e?"6*7** 


oY = 


Similarly, from (43) and the definition of A(v), we have 


0 
J an(o)Wus(0) = -Ha(0), (52) 
-1 
and also from (39), we have 


a) (53) 


Fy1(0) = —g11q(0) — gi q(0) + 7° ( —aRe{M} 


Now we substitute (49) and (53) into (52), and we get 


ef i 7911 , MOLT cats 
(Pt) Bao + Ba) +0) 
fe O re ere 
= (8S Bhar + BA gH" + C2) 


-nareaner (san). 


Thus, we have 


(Fa) a= (att) 
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which gives 


et OY —Re{M} 
2= 0?) ~ \ 4((m— n)Re{M} — 2m) } 
By the above computations, we obtain coefficients of (36) and we can now 
determine the following quantities: 


i 2, 
c1(0) = 3e* (surg 2\ gui \goa| ) = 


3 
Re(c,(0)) 


M2 = ~~ aa(t*) 
Re( A) 


Bg = 2Re(ci (0), 


Ty = _2ileu(0) + potm( Sy). (54) 
&* 

These values describe the bifurcating periodic solutions for the critical value 
T =7*. The direction of the Hopf bifurcation is determined by fz, that 
is, when jig is positive (negative), we have a supercritical (subcritical) Hopf 
bifurcation. For tr > 7* (7 < 7*), there exist bifurcating periodic solutions, 
and {2 determines the stability of the bifurcating periodic solutions: if bg < 
0(682 > 0), then bifurcating periodic solutions in the center manifold are 
stable (unstable). Period of the bifurcating periodic solutions is determined 
by T, the period increases (decrease) when T2 > 0(T) < 0). Thus (20) and 
(54) imply that the results of the direction the Hopf bifurcations satisfy. 


5 Numerical Simulation 


In this section, we study numerical results of system (4) at different values 
of r. By choosing r = 21, w = 21, 6=1, 8 = 3, a = 15 in system (4), we 
have 


d 

a = 21x(1 - sea — 24 — 3xy, 
d 

77 = 3x(t — r)y(t — rT) — 15y. 


For r = 0, Ey = (5,3.75) is an interior unique equilibrium point. These 
values show that (19) has a positive root z = 129.6297769, and then Theorem 
3 shows that 7* = 0.2648100601. From the formulas in Section 4, we get 
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M = —0.3333333334 — 0.7241971227 4, 
N = —0.1467688360 + 1.245730037 i, 
E = —0.3353048563 — 2.701455914 3, 
go = 0.07001726116 + 0.0679380362 1, 
gu = 0.02285330684 — 0.1841225970 1, 
goz = 0.03920150800 — 0.9478826122 i, 


CM = 0.4358591088 + 0.54402174504, 
C? = —0.6386367604 — 0.1582198462 i, 
WY (0) = —0.2638607674 + 0.6954569963:, 
Ww) (0) = —0.3302243769 — 0.5360810164:, 
WY (-1) = 0.1400340404 — 0.2463831080., 
W5o (-1) 
) 


Lip) — 
" = 
2)ig) — 
20 — 
5 (—1) = 1.192302058 — 0.4408369488 i, 


CY = —0.3333333334, 


CY) = —0.4444444444, 
WD (0) = —1.004521703 + 0.03, 
WwW?) (0) = —0.1699352452 + 0.03, 
WP (-1) = —0.9494599988 + 0.04, 


W) (1) = 0.0740685228 + 0.03, 
go. = —0.5523573046 + 1.04184941. 


By replacing the above values in (54), we obtain c,(0) = —0.2658449039 + 
0.46814903207, which implies 2 = 49.95991907 > 0, Bg = —0.5316898078 < 
0 and T> = 0.04309591200 > 0. The fact T> > 0 shows that the period of 
the bifurcating periodic solution increases, and finally, the periodic solution 
disappears. Now, by substituting the values of r = 21, w = 21, d6=1, 8 

3, a= 15 and 7 = 0.2648100601 in system (5) and (6), we can simulate their 
dynamics by changing ¢, numerically. 


In Figure 1, we see that for 7 = 0.05 < r* and ¢ = 100, F = (20,0) in 
(6) and £ = (5,3.75) in (4) and (5) are stable focus. Figure 2 shows that 
for T = 0.2648100601 and ¢ = 0.001, (5) and (6) in F = (20,0) have stable 
focus, but (4) has a periodic solution in F = (5,3.75). In Figure 3, with a 
slight increase in the amount of ¢, no change in the dynamics of (5) and (6) is 
observed. In Figure 4, with increasing ¢ = 10, the stable focus of E = (20,0) 
in (6) remains unchanged. In contrast, the stable focus of (5) changes from 
E = (20,0) to E = (12,3.75), which is a sign of the effectiveness of viral 
therapy and reduction of tumor cells levels. In Figure 5, for ¢ = 100, around 
E = (5,3.75), a family of stable periodic solutions appears in (5) with period 
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less than periodic solutions around F = (5,3.75) in (4), which indicates the 
coexistence between uninfected and infected cells, but H = (20,0) is a stable 
focus in (6). In Figures 6 and 7, we see that with increasing ¢ = 1000 and 
€ = 1000000, a family of stable periodic solutions with the same period will 
appear around £ = (5,3.75), which shows at high values of ¢, (4) and (5) 
have the same dynamic. While in (6), no change in equilibrium point and 
its dynamic is observed and indicate the ineffectiveness of viral therapy in 
reducing the number of infected tumor cells. In Figure 8, for « = 1000 
and 7 = 1, in (4) and (5), the periodic solution disappears. We have used 
Mathematica and Maple software for numerical simulating. 


Gix.y)=8x 


x(t) — er Bil+e)x 
ee (x+ ¥+e) 


Gxy)= 


ixy 3 | 


10 20 30 40 50 


y(t) 
8 


> A 


i) 


t 
10 20 30 40 50 


Figure 1: For « = 100 and r* = 0.05, E = (20,0) in (6), and EF = (5,3.75) in (4) and 
(5) are stable focus. 
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Figure 2: For « = 0.001 and r* = 0.2648100601 E = (20,0) in (5), and (6) is a stable 
focus and in (4) EF = (5,3.75) enclosed by a periodic solution. 
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Figure 3: For ¢ = 1 and r* = 0.2648100601 F = (20,0) in (5) and (6) is a stable focus 


and EF = (5,3.75) in (4) enclosed by a periodic solution. 
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Figure 4: For « = 10 and 7* = 0.2648100601 E = (20,0) in (6) is a stable focus, 
E = (5,3.75) in (4) enclosed by a periodic solution and E = (12,3.75) in (5) is a stable 
focus. 
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Figure 5: For « = 100 and r* = 0.2648100601 E = (20,0) in (6) is a stable focus, 


E = (5,3.75) in (4) and in (5) enclosed by a periodic solution with different period of 
the bifurcating. 
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Figure 6: For « = 1000 and 7* = 0.2648100601 E = (20,0) in (6) is a stable focus, 
E = (5,3.75) in (4) and in (5) enclosed by a periodic solution with almost the same 
period. 


Hopf bifurcation analysis in a delayed model of tumor therapy ... 


Gix.y)=hx 


Gix.y)= 
(x+ ¥+8) 


Gx.y)= 


— 
i=) 


10 20 30 40 £50 


Wy 


t 
10 20 30 40 50 


B(l+e)x 


191 


Figure 7: For ¢ = 1000000 and 7* = 0.2648100601 E = (20,0) in (6) is a stable focus, 


E = (5,3.75) in (4) and (5) enclosed by the same periodic solution. 
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Figure 8: For « = 1000 and 7 = 1 E = (20,0) in (6) is a stable focus, and in (4) and (5) 
the periodic solution disappears. 


6 Discussion 


Three nonlinear mathematical models were expressed with delayed differen- 
tial equations. Two rates of infection, “rapid spread of the virus,” and one 
rate of infection, “slow spread of the virus,” are proposed. Due to the sim- 
plicity of calculating the mathematical model with a linear infection rate, the 
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equilibrium points stability was studied without time delay. By introducing 
a time delay 7 to the second equation in (4), we obtained a periodic solu- 
tion and (4) undergoes a Hopf bifurcation. The existence of periodic solutions 
showed a coexistence between uninfected tumor cells and infected tumor cells, 
leading to a dynamical balance between in growth of uninfected cancer cells 
and infected ones. Now by placing the appropriate value 7 obtained from 
(4) and different values of ¢ in (5) and (6), we simulated the existence of a 
periodic solution, numerically. In the simulation obtained from these three 
models, we found that in (4) and (5) with “rapid virus spread” rate, a peri- 
odic solution can be observed. From a biological point of view, it was stated 
that viral therapy in models with a “rapid virus spread” rate can create a 
coexistence interaction between uninfected and infected tumor cells. While 
simulation (6) showed that no limit cycles are observed for these values, viral 
treatment may not reduce tumor cells. By increasing the time delay, in (4) 
and (5), the period of the created periodic solutions increased until oscilla- 
tions between two populations vanish. Hence, the growth of tumor cancer 
cells became stable. 
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